Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data

We propose a robust interpolation scheme for aliased regularly sampled seismic data that uses the curvelet transform. In a first pass, the curvelet transform is used to compute the curvelet coefficients of the aliased seismic data. The aforementioned coefficients are divided into two groups of scales: alias-free and alias-contaminated scales. The alias-free curvelet coefficients are upscaled to estimate a mask function that is used to constrain the inversion of the alias-contaminated scale coefficients. The mask function is incorporated into the inversion via a minimum norm least squares algorithm that determines the curvelet coefficients of the desired alias free data. Once the aliasfree coefficients are determined, the curvelet synthesis operator is used to reconstruct seismograms at new spatial positions. The proposed method can be used to reconstruct both regularly and irregularly sampled seismic data. We believe that our exposition leads to a clear unifying thread between f-x Spitz (1991) and f-k Gulunay (2003) beyond-alias interpolation methods and curvelet reconstruction. Likewise in f-x and f-k interpolation, we stress the necessity of examining seismic data at different scales (frequency bands) in order to come up with viable and robust interpolation schemes. Synthetic and real data examples are used to illustrate the performance of the proposed curvelet interpolation method. INTRODUCTION Interpolation and reconstruction of seismic data has become an important topic for the seismic processing community. It is often the case that logistic and economic constraints dictate the spatial sampling of seismic surveys. Wave fields are continuous, in other words, seismic energy reaches the surface of the earth everywhere in our area of study. The process of acquisition records a finite number of spatial samples of the continuous wave field generated by a finite number of sources. The latter leads to a regular or irregular distribution of sources and receivers. Many important techniques for removing coherent noise and imaging the earth interior have stringent sampling requirements which are often not met in real surveys. In order to avoid information losses, the data should be sampled according to the Nyquist criterion (Vermeer, 1990). When this criterion is not honored, reconstruction can be used to recover the data to a denser distribution of sources and receivers and mimic a properly sampled survey (Liu, 2004). The final result of the reconstruction stage could have a significant impact on subsequent seismic processing steps such as noise removal (Soubaras, Curvelet Interpolation 2 2004), AVO analysis (Sacchi and Liu, 2005; Hunt et al., 2008), and imaging (Liu and Sacchi, 2004). Methods for seismic wave field reconstruction can be classified into two categories: waveequation based methods and signal processing methods. Wave-equation methods utilize the physics of wave propagation to reconstruct seismic volumes. In general, the idea can be summarized as follows. An operator is used to map seismic wave fields to a physical domain. Then, the modeled physical domain is transformed back to data space to obtain the data we would have acquired with an ideal experiment. It is basically a regression approach where the regressors are built based on wave equation principles (in general, approximations to kinematic ray theoretical solutions of the wave equation). The methods proposed by Ronen (1987), Bagaini and Spagnolini (1999), Stolt (2002), Trad (2003), Fomel (2003), Malcolm et al. (2005), Clapp (2006) and Leggott et al. (2007) fall under this category. These methods require the knowledge of some sort of velocity distribution in the earth’s interior (migration velocities, rms velocities, stacking velocities). While reconstruction methods based on wave equation principles are very important, this paper will not investigate this category of reconstruction algorithms. Seismic data reconstruction via signal processing approaches is an ongoing research topic in exploration seismology. During the last decade, important advances have been made in this area. Nowadays, signal processing reconstruction algorithms based on Fourier synthesis operators can cope with multidimensional sampling as demonstrated by several authors (Duijndam et al., 1999; Liu et al., 2004; Zwartjes and Gisolf, 2006; Schonewille et al., 2009). These methods are based on classical signal processing principles, they do not require information about the subsurface and, in addition, are quite robust in situations were the optimality condition under which they were designed are not completely satisfied (Trad, 2009). Signal processing methods for seismic data reconstruction often rely on transforming the data to other domains. The most commonly used transformations are the Fourier transform (Sacchi and Ulrych, 1996; Sacchi et al., 1998; Duijndam et al., 1999; Liu et al., 2004; Xu et al., 2005; Zwartjes and Gisolf, 2006), the Radon transform (Darche, 1990; Verschuur and Kabir, 1995; Trad et al., 2002), the local Radon transform (Sacchi et al., 2004; Wang et al., 2009) and the curvelet transform (Hennenfent and Herrmann, 2008; Herrmann and Hennenfent, 2008). Another group of signal processing interpolation methods rely on prediction error filtering techniques (Wiggins and Miller, 1972). Spitz (1991) and Porsani (1999) introduced beyond-alias seismic trace interpolation methods using prediction filters. These methods operate in the frequency-space (f-x ) domain. In both cases, low frequency data components in a regular spatial grid are used to estimate the prediction filters needed to interpolate high frequency data components. An equivalent interpolation method in the frequency-wavenumber (f-k) domain was introduced by Gulunay (2003) and often referred as f-k interpolation. The main contribution of this paper is the introduction of a strategy that utilizes the curvelet transform to interpolate regularly sampled aliased seismic data. It is important to stress that the curvelet transform has been used by Hennenfent and Herrmann (2007), Hennenfent and Herrmann (2008) and Herrmann and Hennenfent (2008) to interpolate seismic data. In their articles, they reported the difficulty of interpolating regularly sampled aliased data with the curvelet transform and therefore, proposed random sampling strategies to circumvent the aliasing problem. This paper, however, proposes a new methodology which successfully eliminates the requirement of randomization to avoid aliasing. We create 3 a mask function from the alias-free curvelet scales (low frequencies) to constrain the interpolation of alias-contaminated scales (high frequencies). The proposed method is an attempt to utilize early principles of f-x and f-k domain beyond alias interpolation methods (Spitz, 1991; Gulunay, 2003) in the curvelet domain. In summary, by carefully understanding wellestablished beyond alias interpolation methods (Spitz, 1991; Gulunay, 2003) in conjunction with novel signal processing tools like the curvelet transform (Candes et al., 2005), we were able to develop an algorithm capable of reconstructing aliased regularly sampled data. In addition, the proposed algorithm can deal with the simpler problem of regularization of randomly sampled data where the alias footprint is annihilated by random sampling. In view of the fact that the curvelet transform is a local transformation, the proposed algorithm can easily cope with strong variations of dips. This is a clear advantage with respect to reconstruction via Fourier bases (Duijndam et al., 1999; Liu, 2004; Liu et al., 2004; Zwartjes and Gisolf, 2006; Trad, 2009) and prediction filtering techniques (Spitz, 1991; Naghizadeh and Sacchi, 2009b) that assume data composed of a superposition of a few plane waves. The latter requires spatial and temporal windowing to satisfy the plane wave assumption. One needs to stress, however, that efficient multidimensional algorithms based on Fourier bases and prediction filtering are already part of current industrial methodologies to reconstruct data that depend on 3 and 4 spatial dimensions (Schonewille et al., 2003; Abma and Kabir, 2006; Schonewille et al., 2009; Trad, 2009). Embarking on multidimensional reconstruction by means of curvelet basis functions could be extremely difficult as, to our knowledge, multidimensional curvelet transforms are limited to 3 dimensions (Ying et al., 2005). Furthermore, curvelet transform reconstruction algorithms appear, at present time, not to be amenable of a fast implementation such as the one encountered in 4D spatial Fourier reconstruction algorithms that solely rely on the Fast Fourier Transform (Liu, 2004; Trad, 2009). We stress that we are not attempting to contrast curvelet based reconstruction with well-tested efficient multidimensional Fourier reconstruction algorithms. Our exposition mainly aims at presenting an interesting strategy for reconstructing aliased data using the curvelet transform. It is clear that an important amount of work is needed to implement curvelet regularization strategies that can be used to regularize 4D spatial data. THEORY The curvelet transform The curvelet transform is a local and directional decomposition of an image (data) into harmonic scales (Candes and Donoho, 2004). The curvelet transform aims to find the contribution from each point of data in the t-x to isolated directional windows in the f-k domain. If we assume that m(t, x) represents seismic data in the t-x domain, we define a set of basis functions as φ(s, θ, t, x), where s indicates scale (increasing from coarsest to finest), θ, is angle or dip and t0, x0 are the t-x location parameters. These basis functions, known as curvelet basis functions, are used to decompose the original data in local components of various scales and dip. Curvelets can be considered as wavelets with the additional important property of directionality (dip). The continuous curvelet transform can be represented as the inner product of the data m(t, x) and curvelet basis functions c(s, θ, t0, x0) = C[m] = ∫ t ∫ x m(t, x)φ(s, θ, t0 − t, x0 − x)dtdx . (1)